% The code is based on the following paper:
%   Chen, T. Y., Lin, Y. L., and Tzeng, L. Y., forthcoming.
%   Estimating probability weighting functions through option pricing bounds. Review of Asset Pricing Studies.
%
% Copyright: Tzu-Ying Chen, Yo-Lan Lin, Larry Y. Tzeng
% Date: January 30, 2024

clear

%% Setting
RootName_DIST = ['EGARCH'];
% RootName_DIST = ['COND_ATMIV'];
% RootName_DIST = ['COND_EGARCHV'];

%% Load Data
FileName = ['AllTradingDate19962021_Monthly.mat'];
load(FileName, ...
     'Target_AllDate');             
clear FileName

Date_Monthly = Target_AllDate;
clear Target_AllDate 

%% Calculate Coverage Rate for Each Parameters
% Transaction Cost
TC = 0.005;

EST_RDEU = nan(size(Date_Monthly, 1), 1);                                  
for d = 1:size(Date_Monthly, 1)
    Target_Date = Date_Monthly(d, 1);
    Target_TTM = Date_Monthly(d, 3);

    % Load Data 
    FileName = ['OP' num2str(Target_Date) '_TTM' num2str(Target_TTM) '.mat'];
    load(FileName, ...
         'Data', ...
         'Index_Date', 'Index_TTM', 'Index_CPFlag', 'Index_K', 'Index_OP_Bid', 'Index_OP_Ask', ...
         'Index_RF', 'Index_DY', 'Index_IS', 'Index_IS_ADJ', 'Index_IV_Bid', 'Index_IV_Ask');
    clear FileName 

    Index_Drop = find(sum(isnan(Data(:, [Index_IV_Bid Index_IV_Ask])), 2) > 0);
    Data(Index_Drop, :) = [];                                              
    AllData = Data;
    clear Index_Drop Data

    K = AllData(:, Index_K);
    S0 = AllData(:, Index_IS);
    Index_KS = size(AllData, 2) + 1;
    AllData(:, Index_KS) = K ./ S0;                                        
    AllData = sortrows(AllData, Index_KS);                                  
    clear K S0

    % Gross Return Distributions
    if strcmp(RootName_DIST, 'EGARCH')==1
        FileName = ['DG_GrossRET_Monthly_' num2str(Target_Date)];
        load(FileName, ...
             'RET_Monthly');
        clear FileName

        [AllRET, ~, Index_AllRET] = unique(RET_Monthly); 
        AllRET_PROP = accumarray(Index_AllRET, 1, size(AllRET), @sum, nan) / length(RET_Monthly);
        State_Monthly = AllRET';
        State_PROB = AllRET_PROP';
        clear Index_AllRET RET_Monthly AllRET AllRET_PROP
    else
    end
    
    if strcmp(RootName_DIST, 'COND_ATMIV')==1
        FileName = ['CJP_ATMIV_GrossRET_Monthly_' num2str(Target_Date)];
        load(FileName, ...
             'State_Monthly', 'State_PROB');                     
        clear FileName        
    else
    end    
        
    if strcmp(RootName_DIST, 'COND_EGARCHV')==1
        FileName = ['CJP_EGARCHV_GrossRET_Monthly_' num2str(Target_Date)];
        load(FileName, ...
             'State_Monthly', 'State_PROB');                     
        clear FileName        
    else
    end      

    % Risk Preference Parameters Estimation
    K = AllData(:, Index_K);
    KS = AllData(:, Index_KS);
    OP_Bid = AllData(:, Index_OP_Bid);
    OP_Ask = AllData(:, Index_OP_Ask);
    CPFlag = AllData(:, Index_CPFlag);
    IV_Bid = AllData(:, Index_IV_Bid);
    IV_Ask = AllData(:, Index_IV_Ask);
    S0 = AllData(:, Index_IS);
    S0_ADJ = AllData(:, Index_IS_ADJ);
    TTM = AllData(:, Index_TTM) / 365;                                     
    RF = AllData(:, Index_RF);                                             
    DY = AllData(:, Index_DY);                                             

    % Expected Utility                  
    OP_LB_EU = RDEU_OP_LB(State_Monthly, State_PROB, ...
                          K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                          1, 1);                
    OP_UB_EU = RDEU_OP_UB(State_Monthly, State_PROB, ...
                          K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                          1, 1);   

    % Rank-Dependent Expected Utility
    if (sum(max(OP_LB_EU, 0) > OP_Ask') > 0) | (sum(max(OP_UB_EU, 0) < OP_Bid') > 0)
        dx = 0.01;
        x1 = 0:dx:3;
        x1(x1==0) = [];                                                    
        x2 = 0:dx:3; 
        x2(x2==0) = [];                                                    
        x1 = repmat(x1, length(x2), 1);                                    
        x2 = repmat(x2', 1, size(x1, 2));                                  
        x = [x1(:) x2(:)];                                                 
        clear dx x1 x2
        
        Record = nan(size(x, 1), size(x, 2) + 3);                          
        Record_OP_LB = nan(size(x, 1), length(K));                         
        Record_OP_UB = nan(size(x, 1), length(K));                         
        for i = 1:size(x, 1)
            OP_LB_RDEU = RDEU_OP_LB(State_Monthly, State_PROB, ...
                                    K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                                    x(i, 1), x(i, 2));   
            OP_LB = min(OP_LB_RDEU, OP_LB_EU);    

            OP_UB_RDEU = RDEU_OP_UB(State_Monthly, State_PROB, ...
                                    K, CPFlag, TTM, S0, S0_ADJ, RF, DY, TC, ...
                                    x(i, 1), x(i, 2));   
            OP_UB = max(OP_UB_RDEU, OP_UB_EU);  

            CR = sum(OP_LB <= OP_Ask', 2);
            CR_LB = CR / length(K);                                        
            clear CR
    
            CR = sum(OP_UB >= OP_Bid', 2);
            CR_UB = CR / length(K);                                        
            clear CR
    
            CR = sum(OP_UB >= OP_Bid', 2) + ...
                 sum(OP_LB <= OP_Ask', 2);
            CR_LBUB = CR / (2 * length(K));                                
            clear CR
    
            Record(i, :) = [x(i, :) CR_LB CR_UB CR_LBUB];
            Record_OP_LB(i, :) = OP_LB;
            Record_OP_UB(i, :) = OP_UB;
            clear OP_LB_RDEU OP_UB_RDEU OP_LB OP_UB CR_LB CR_UB CR_LBUB
        end
        clear x
        
        EST_RDEU(d, :) = 1;                                                
    else
        x = [1 1];
        CR_LB = 1;                                                         
        CR_UB = 1;                                                         
        CR_LBUB = 1;                                                       
        Record = [x CR_LB CR_UB CR_LBUB];
        Record_OP_LB = OP_LB_EU;
        Record_OP_UB = OP_UB_EU;        
        clear x CR_LB CR_UB CR_LBUB
        
        EST_RDEU(d, :) = 0;                                                
    end
    
    % Output
    FileName = ['Summary_OP_Grid_' num2str(Target_Date)];
    save(FileName, ...
         'Record', 'Record_OP_LB', 'Record_OP_UB', ...
         'OP_LB_EU', 'OP_UB_EU', ...
         'AllData', ...
         'OP_Bid', 'OP_Ask', 'IV_Bid', 'IV_Ask', ...
         'K', 'KS', 'CPFlag', 'S0', 'S0_ADJ', 'TTM', 'RF', 'DY', ...
         'TC', ...
         'State_Monthly', 'State_PROB');
    clear FileName    
    
    clear Target_Date Target_TTM AllData State_Monthly State_PROB
    clear Index_Date Index_TTM Index_CPFlag Index_K Index_OP_Bid Index_OP_Ask
    clear Index_RF Index_DY Index_IS Index_IS_ADJ Index_IV_Bid Index_IV_Ask Index_KS 
    clear K KS CPFlag OP_Bid OP_Ask IV_Bid IV_Ask S0 S0_ADJ TTM RF DY
    clear OP_LB_EU OP_UB_EU Record Record_OP_LB Record_OP_UB
end
clear TC
clear d i